subroutine faraday
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: dummy
!
!
      bxv(:ibp2*jbp2*kbp2) = 0.
      byv(:ibp2*jbp2*kbp2) = 0.
      bzv(:ibp2*jbp2*kbp2) = 0.
!
!     ******************************************************
!
!     solve Faraday's law to advance the magnetic field (cell centered)
!
         call curlc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, curlx, curly, curlz)
!
!
!	apply bc on B cell centered
!
!	Note that the last entry is set to 0 to impose zero tangential value
!       in the ghost and boundary centres, 1 to apply neumann there. Zero value needs
!	to be applied to the curl E but 1d0 needs to be applied to full B
!


      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         curlx, curly, curlz, periodicx,0d0)
!
         dtc = clite*dt
!
!         do n = 1, ncells
!
!            ijk = ijkcell(n)
!
         do i=1,ibp2
         do j=1,jbp2
         do k=1,kbp2

            ijk= 1+ (i-1)*iwid+(j-1)*jwid+(k-1)*kwid

            bxn(ijk) = bxn(ijk) - curlx(ijk)*dtc
            byn(ijk) = byn(ijk) - curly(ijk)*dtc
            bzn(ijk) = bzn(ijk) - curlz(ijk)*dtc
!
         end do
         end do
         end do
!
!     ****************************************************
!
!	apply bc on B cell centered
!
      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx,1d0)
!
      call vtxb (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, vol, bxn, byn, &
         bzn, bxv, byv, bzv)
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, wate&
         (1,22))
!
      call bcperv (ibp2, jbp2, kbp2, iwid, jwid, kwid, sdlt, cdlt, bxv, byv, &
         bzv, periodicx)

      call boundary_b_vertex(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxv, byv, bzv, periodicx,1d0)
!
      call current

!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)
!
      return
!
!	The stuff below is an ancient relic.
!
!      call bdrift (ncells, ijkcell, nsp, itdim, iwid, jwid, kwid, qdnv, qom, &
!         bxn, byn, bzn, dt, clite, wate(1,1), wate(1,2), wate(1,3))
!
!     impose toroidal boundary conditions on "alfab"
!
!      dummy = 0.0
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, wate(1,1), wate(1,&
!         2), wate(1,3), periodicx)
!
!      call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
!         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
!         , c5z, c6z, c7z, c8z, wate(1,1), wate(1,2), wate(1,3), udrft, vdrft, &
!         wdrft)
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, udrft, vdrft, &
!         wdrft, periodicx)
!
!      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, udrft, vdrft, wdrft)
!
!      write (*, *) 'vinit:renormalize udrift'
!      do n = 1, nvtx
!         ijk = ijkvtx(n)
!
!         rvvol = 1./vvol(ijk)
!
!         udrft(ijk) = udrft(ijk)*rvvol
!         vdrft(ijk) = vdrft(ijk)*rvvol
!         wdrft(ijk) = wdrft(ijk)*rvvol
!
!      end do
!
      end subroutine faraday

subroutine faraday_new
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      use boundary
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: dummy, dtcth
    real(double) :: density,temperature,flvx,flvy,flvz,bx0_ini, &
                    by0_ini,bz0_ini, bx1_ini,by1_ini,bz1_ini
!
!     ******************************************************
!
!
      bxv(:ibp2*jbp2*kbp2) = 0.
      byv(:ibp2*jbp2*kbp2) = 0.
      bzv(:ibp2*jbp2*kbp2) = 0.
!
      dtcth = clite*dt*cntr
!
      call current
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)

      jx0old(:ibp2*jbp2*kbp2) = jx0(:ibp2*jbp2*kbp2)
      jy0old(:ibp2*jbp2*kbp2) = jy0(:ibp2*jbp2*kbp2)
      jz0old(:ibp2*jbp2*kbp2) = jz0(:ibp2*jbp2*kbp2)
!
!
!     solve Faraday's law to advance the magnetic field (cell centered)
!

         call curlc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, J12x, J12y, J12z, curlx, curly, curlz)
!
!     impose toroidal boundary conditions on curl(J)
!

      dummy = 0.0D0

      call torusbc (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, curlx, &
         curly, curlz, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, curlx, &
         curly, curlz)
!

         call curlc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0, wate(1,1), wate(1,2), wate(1,3))
!
!     impose toroidal boundary conditions on curl(E)
!

      dummy = 0.0D0

      call torusbc (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, &
          wate(1,1), wate(1,2), wate(1,3), periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, wate(1,1), wate(1,2), wate(1,3))
!
!
         curlx(:ibp2*jbp2*kbp2)=-curlx(:ibp2*jbp2*kbp2)*(dtcth)**2 &
                                +wate(:ibp2*jbp2*kbp2,1)*dtcth-bxn(:ibp2*jbp2*kbp2)
         curly(:ibp2*jbp2*kbp2)=-curly(:ibp2*jbp2*kbp2)*(dtcth)**2 &
                                +wate(:ibp2*jbp2*kbp2,2)*dtcth-byn(:ibp2*jbp2*kbp2)
         curlz(:ibp2*jbp2*kbp2)=-curlz(:ibp2*jbp2*kbp2)*(dtcth)**2 &
                                +wate(:ibp2*jbp2*kbp2,3)*dtcth-bzn(:ibp2*jbp2*kbp2)
!
!
write(*,*)'Norm Bx (voor)=',sum(bxn(:ibp2*jbp2*kbp2)**2),sum(curlx(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm By (voor)=',sum(byn(:ibp2*jbp2*kbp2)**2),sum(curly(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Bz (voor)=',sum(bzn(:ibp2*jbp2*kbp2)**2),sum(curlz(:ibp2*jbp2*kbp2)**2)

!
call boundary_state(ibar/2,jbar/2,2,density,temperature,flvx,flvy,flvz,bx0_ini,by0_ini,bz0_ini)
call boundary_state(ibar/2,jbar/2,kbar+1,density,temperature,flvx,flvy,flvz,bx1_ini,by1_ini,bz1_ini)

         rdt=-1.d0
         scalsq=-dtcth**2
         wk = 0.
         phibc = 1.d0
         phibce = bx0_ini
         phibcf = bx1_ini
         write(*,*)'bx0',bx0_ini,bx1_ini
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlx, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), bxn, wate(1,7), wkl, wkr, wkf, wke, periodicx)

         ijk = 1 + (ibar - 1)*iwid + (jbar - 1)*jwid + kbp1*kwid
         write(*,*)'*******',bxn(ijk),bxn(ijk-kwid)
                  ijk = 1 + (ibar - 1)*iwid + (jbar - 1)*jwid
         write(*,*)'*******',bxn(ijk),bxn(ijk+kwid)
!
         phibc=1d0
         phibce = by0_ini
         phibcf = by1_ini
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curly, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), byn, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         phibc=1d0
         phibce = bz0_ini
         phibcf = bz1_ini
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlz, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), bzn, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
write(*,*)'Norm Bx (na)=',sum(bxn(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm By (na)=',sum(byn(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Bz (na)=',sum(bzn(:ibp2*jbp2*kbp2)**2)

      wk = 1.
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, bxn, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, byn, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, bzn, phibc, wk, a11, periodicx)
!
!
      call boundary_b_vertex(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxv, byv, bzv, periodicx,1d0)
!
      call current
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)
!
         curlx(:ibp2*jbp2*kbp2)=-(-jx0old(:ibp2*jbp2*kbp2)+jx0(:ibp2*jbp2*kbp2))*clite/fourpi
         curly(:ibp2*jbp2*kbp2)=-(-jy0old(:ibp2*jbp2*kbp2)+jy0(:ibp2*jbp2*kbp2))*clite/fourpi
         curlz(:ibp2*jbp2*kbp2)=-(-jz0old(:ibp2*jbp2*kbp2)+jz0(:ibp2*jbp2*kbp2))*clite/fourpi
!
write(*,*)'Norm Ex (voor)=',sum(ex(:ibp2*jbp2*kbp2)**2),sum(curlx(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ey (voor)=',sum(ey(:ibp2*jbp2*kbp2)**2),sum(curly(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ez (voor)=',sum(ez(:ibp2*jbp2*kbp2)**2),sum(curlz(:ibp2*jbp2*kbp2)**2)

!
goto 666
!
         rdt=0d0
         scalsq=dtcth
         wk = 0.
         phibc = 0.d0
         wke=1d0
         wkf=1d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlx, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ex, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         phibc=0d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curly, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ey, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         phibc=0d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlz, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ez, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         wke=0d0
         wkf=0d0
!
write(*,*)'Norm Ex (na)=',sum(ex(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ey (na)=',sum(ey(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ez (na)=',sum(ez(:ibp2*jbp2*kbp2)**2)
!

666      return
!
!	The stuff below is an ancient relic.
!
!      call bdrift (ncells, ijkcell, nsp, itdim, iwid, jwid, kwid, qdnv, qom, &
!         bxn, byn, bzn, dt, clite, wate(1,1), wate(1,2), wate(1,3))
!
!     impose toroidal boundary conditions on "alfab"
!
!      dummy = 0.0
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, wate(1,1), wate(1,&
!         2), wate(1,3), periodicx)
!
!      call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
!         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
!         , c5z, c6z, c7z, c8z, wate(1,1), wate(1,2), wate(1,3), udrft, vdrft, &
!         wdrft)
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, udrft, vdrft, &
!         wdrft, periodicx)
!
!      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, udrft, vdrft, wdrft)
!
!      write (*, *) 'vinit:renormalize udrift'
!      do n = 1, nvtx
!         ijk = ijkvtx(n)
!
!         rvvol = 1./vvol(ijk)
!
!         udrft(ijk) = udrft(ijk)*rvvol
!         vdrft(ijk) = vdrft(ijk)*rvvol
!         wdrft(ijk) = wdrft(ijk)*rvvol
!
!      end do
!
      end subroutine faraday_new


subroutine faraday_new_new
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      use boundary
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: dummy, dtcth
    real(double) :: density,temperature,flvx,flvy,flvz,bx0_ini, &
                    by0_ini,bz0_ini, bx1_ini,by1_ini,bz1_ini
!
!     ******************************************************
!
!
      bxv(:ibp2*jbp2*kbp2) = 0.
      byv(:ibp2*jbp2*kbp2) = 0.
      bzv(:ibp2*jbp2*kbp2) = 0.
!
      dtcth = clite*dt*cntr
!
      call current
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)

      jx0old(:ibp2*jbp2*kbp2) = jx0(:ibp2*jbp2*kbp2)
      jy0old(:ibp2*jbp2*kbp2) = jy0(:ibp2*jbp2*kbp2)
      jz0old(:ibp2*jbp2*kbp2) = jz0(:ibp2*jbp2*kbp2)
!
!
      bxv(:ibp2*jbp2*kbp2) = 0.
      byv(:ibp2*jbp2*kbp2) = 0.
      bzv(:ibp2*jbp2*kbp2) = 0.
!
!     ******************************************************
!
!     solve Faraday's law to advance the magnetic field (cell centered)
!
         call curlc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, curlx, curly, curlz)
!
!
!	apply bc on B cell centered
!
!	Note that the last entry is set to 0 to impose zero tangential value
!       in the ghost and boundary centres, 1 to apply neumann there. Zero value needs
!	to be applied to the curl E but 1d0 needs to be applied to full B
!


      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         curlx, curly, curlz, periodicx,0d0)
!
         dtc = clite*dt
!
!         do n = 1, ncells
!
!            ijk = ijkcell(n)
!
         do i=1,ibp2
         do j=1,jbp2
         do k=1,kbp2

            ijk= 1+ (i-1)*iwid+(j-1)*jwid+(k-1)*kwid

            bxn(ijk) = bxn(ijk) - curlx(ijk)*dtc
            byn(ijk) = byn(ijk) - curly(ijk)*dtc
            bzn(ijk) = bzn(ijk) - curlz(ijk)*dtc
!
         end do
         end do
         end do
!
!     ****************************************************
!
!	apply bc on B cell centered
!
      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx,1d0)
!
      call vtxb (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, vol, bxn, byn, &
         bzn, bxv, byv, bzv)
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, wate&
         (1,22))
!
      call bcperv (ibp2, jbp2, kbp2, iwid, jwid, kwid, sdlt, cdlt, bxv, byv, &
         bzv, periodicx)

      call boundary_b_vertex(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxv, byv, bzv, periodicx,1d0)
!
      call current

!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)
!
!
      call current
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         jx0, jy0, jz0, periodicx)
!
         curlx(:ibp2*jbp2*kbp2)=-(-jx0old(:ibp2*jbp2*kbp2)+jx0(:ibp2*jbp2*kbp2))*clite/fourpi
         curly(:ibp2*jbp2*kbp2)=-(-jy0old(:ibp2*jbp2*kbp2)+jy0(:ibp2*jbp2*kbp2))*clite/fourpi
         curlz(:ibp2*jbp2*kbp2)=-(-jz0old(:ibp2*jbp2*kbp2)+jz0(:ibp2*jbp2*kbp2))*clite/fourpi
!
write(*,*)'Norm Ex (voor)=',sum(ex(:ibp2*jbp2*kbp2)**2),sum(curlx(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ey (voor)=',sum(ey(:ibp2*jbp2*kbp2)**2),sum(curly(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ez (voor)=',sum(ez(:ibp2*jbp2*kbp2)**2),sum(curlz(:ibp2*jbp2*kbp2)**2)

!
!goto 666
!
         rdt=0d0
         scalsq=dtcth
         wk = 0.
         phibc = 0.d0
         wke=1d0
         wkf=1d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlx, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ex, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         phibc=0d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curly, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ey, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         phibc=0d0
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,8), wate(1,9)&
                  , wate(1,10), wate(1,11), itmax, error, curlz, rdt, &
                  scalsq, wk, phibc, &
                  wate(1,1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                  6), ez, wate(1,7), wkl, wkr, wkf, wke, periodicx)
!
         wke=0d0
         wkf=0d0
!
write(*,*)'Norm Ex (na)=',sum(ex(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ey (na)=',sum(ey(:ibp2*jbp2*kbp2)**2)
write(*,*)'Norm Ez (na)=',sum(ez(:ibp2*jbp2*kbp2)**2)
!

666      return
!
!	The stuff below is an ancient relic.
!
!      call bdrift (ncells, ijkcell, nsp, itdim, iwid, jwid, kwid, qdnv, qom, &
!         bxn, byn, bzn, dt, clite, wate(1,1), wate(1,2), wate(1,3))
!
!     impose toroidal boundary conditions on "alfab"
!
!      dummy = 0.0
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, wate(1,1), wate(1,&
!         2), wate(1,3), periodicx)
!
!      call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
!         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
!         , c5z, c6z, c7z, c8z, wate(1,1), wate(1,2), wate(1,3), udrft, vdrft, &
!         wdrft)
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, udrft, vdrft, &
!         wdrft, periodicx)
!
!      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, udrft, vdrft, wdrft)
!
!      write (*, *) 'vinit:renormalize udrift'
!      do n = 1, nvtx
!         ijk = ijkvtx(n)
!
!         rvvol = 1./vvol(ijk)
!
!         udrft(ijk) = udrft(ijk)*rvvol
!         vdrft(ijk) = vdrft(ijk)*rvvol
!         wdrft(ijk) = wdrft(ijk)*rvvol
!
!      end do
!
      end subroutine faraday_new_new

subroutine div_cleaning
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
	real(double) :: rcntr, turnon
!
!
!     advance electric field from E(t+cntr*dt) to E(t+dt)
!
      rcntr = 1./cntr
      if (eandm) then
         turnon = 1.
      else
         turnon = 0.
      endif
!
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
!
         ex(ijk) = (ex0(ijk)+(ex(ijk)-ex0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ex(ijk)
         ey(ijk) = (ey0(ijk)+(ey(ijk)-ey0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ey(ijk)
         ez(ijk) = (ez0(ijk)+(ez(ijk)-ez0(ijk))*rcntr)*turnon + (1. - turnon)*&
            ez(ijk)
!
         ex0(ijk) = ex(ijk)
         ey0(ijk) = ey(ijk)
         ez0(ijk) = ez(ijk)
!
      end do

!     calculate the source for Poisson's equation
!     the call computes a correction to E0
!
               call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
	dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
		  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'prima=',dnorm
!
!     with these parameters, poisson solves poisson's equation
!
               dumdt = 0.0d0
               dumscal = -1.d0 !-1 for poisson not fft
               wk = 0.
               phibc = 0.d0
!
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
!
               call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
                  c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
                  c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, phipot, gradxb, &
                  gradyb, gradzb)
!
!
               call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, &
                  gradzb) !....nothing is done with this call
!
!dir$ ivdep
               do n = 1, nvtx
                  ijk = ijkvtx(n)
                  rvvol = 1./vvol(ijk)
                  ex0(ijk) = ex0(ijk) - gradxb(ijk)*rvvol
                  ey0(ijk) = ey0(ijk) - gradyb(ijk)*rvvol
                  ez0(ijk) = ez0(ijk) - gradzb(ijk)*rvvol
               enddo

         call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
			   dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
				  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'dopo=',dnorm



end subroutine div_cleaning

